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Abstract. We report on the development of an implicit multi-D hydrody- 
namic code for stellar evolution. We present two test-cases relevant for the first 
scientific goal of the code: the simulation of convection in pulsating stars. First 
results on a realistic stellar model are also presented. 



1. Introduction 

Time explicit schemes are prone to the Courant-Friedriech-Lewy (CFL) stabil- 
ity condition which constrains the maximum time step one can use in numerical 
calculations. However, in some stellar physics problems, the physical solution 
evolves on time-scales which are much longer than this "numerical" time-scale. 
This often makes calculations too cumbersome to be tractable. It is then neces- 
sary to use an implicit scheme, which does not have such a limitation on the time 
step. This paper describes such an implicit hydrodynamic code. The first scien- 
tific problem to be adressed with the code is the role of convection in pulsating 
stars. This requires an implicit approach since the typical CFL time step (for a 
moderate resolution) is AtcFL ~ 10 s, whereas the convection (eddy turnover) 
time-scale is Tconv ~ 10^~^ s (uconv ~ 0.1 — O.Scsound) and typical pulsation pe- 
riods are of the order of a few days (~ 10^~^ s) with a growth rate of the order 
of ~ 10^ s. 



2. Code description 

We solve numerically the standard hydrodynamic conservation equations with 
radiative transfer. We use the internal energy density equation and we treat 
radiation transport in the diffusion limit. This approximation is suitable for 
optically thick region but becomes inaccurate in the photosphere and optically 
thin regions. Molecular viscosity is negligible in astrophysical flow, but we use 
artificial viscosity when robustness might be an issue, e.g. for the initial relax- 
ation of stellar models or for shocks. The hydrodynamic equations are closed 
with an equation of state and opacities appropriate for the description of stellar 
interior structure. 

We use spherical coordinates and assume azimuthal symmetry, leaving (r, 6) 
as independent coordinates. The spatial discretization is done on a staggered 
grid, using a finite volume approach. Fluxes are computed at cell interfaces 
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Table 1. Run parameters for the entropy bubble test 



Run 


Time 


Resolution 


time step 


CFL 


1 


50 


50^ 


0.04 


1 


2 


50 


502 


0.5 


12.5 


3 


200 


1502 


0.5 


37.5 



with an upwind reconstruction scheme: we use either the donor cell method 
(first order in space) or the Van Leer method (second order in space). For 
the implicit temporal discretization, we use either the backward Euler scheme 
(first order in time) or the Crank-Nicholson scheme (second order in time). An 
implicit discretization results in a non-linear system which is solved by Newton- 
Raphson iterations. This requires the inversion of a linear system involving the 
jacobian matrix of the set of equations. Here the jacobian matrix is computed 
by numerical differencing and we use a direct solver, MUMPS (Amestoy et al. 
2001), which implements an efficient and robust LU decomposition. The variable 
time step is chosen on empirical ground, depending on the amount of changes 
of variables during a time step, as described for example in Dorfi (1997). 



3. Tests cases 

For test purpose, we have also implemented in the code cartesian geometry 
(2D) and the equation of state of a perfect gas with a constant adiabatic index 
7. Basic tests (e.g. Sedov, Sod, Barenblatt) have been successfully done for code 
validation during its first stages of development, so we focus here on two tests 
more relevant to stellar physics problems such as convection or oscillations. 

3.1. Test 1: oscillating entropy bubble 

This test is taken from Dintrans & Brandenburg 2004 (DB in the following) . We 
consider an isothermal stratified atmosphere with constant gravity perturbed by 
an entropy bubble. We compute the oscillations of this bubble around its equilib- 
rium position and we analyze the spectrum of internal gravity and sound waves 
excited by the bubble. The initial setup, parameters and units normalization 
are similar to DB and are therefore not reproduced here. The Brunt- Vaisala 
frequency of the atmosphere is N ~ 0.82. The kinematic viscosity v is set to 
5 X 10"^ in ah runs. The thermal conductivity kr is set to zero, as we are in- 
terested in the adiabatic evolution of the system. We use a cartesian domain 
{x,z) € [—0.5,0.5] X [0,1]. Periodic boundary conditions are used in the hori- 
zontal directions and non penetrative conditions (zero vertical velocity) are used 
at the bottom and top boundaries. All the run parameters are summarized in 
Tabled! 

As in DB, we describe each wave mode by two integers l,n >0 which are 
the number of nodes in the horizontal and vertical direction. In this test we look 
at two types of waves (for the derivation of these properties we refer the reader 
to DB): 1) Vertical sound waves which have frequencies uj = {n + l)7r. These 
waves have I = since they have no horizontal structure. 2) Internal gravity 
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Figure 1. Accoustic modes in the {z,uj) plane for At — 0.04 (left) and 
At = 0.5 (right). 



waves which have frequencies uj < N and are characterized by / 7^ (internal 
gravity modes cannot propagate in the vertical direction). 

These two kinds of waves occupy well separated regions of the frequency 
spectrum. Vertical sound waves have vertical periods P < 2 and will therefore 
be more sensitive to the time step. Internal gravity modes have longer periods 
P > 7.3 and will therefore be less sensitive. This test thus provides a good 
analysis of the influence of the time step on the accuracy of the numerical results. 

We first analyze the vertical sound wave spectrum by using method 2 de- 
scribed in DB, with run 1 &: 2 (c.f. Tabled]). In this method, the wave spectrum 
is represented in the {z, uj) plane. Note that for a given time step At, the Nyquist 
theorem states that no periodic signal with frequency greater than = 7r/At 
can be resolved. Also the temporal Fourier transform has a frequency resolution 
which is Auj = 27r/T, where T is the duration of the simulation. The results 
are shown in fig. [1] for two different time steps. For At = 0.04 we recognize 
the signature of the n = 0, 1, 2 modes (the n = 3 mode is below the range of 
the isocontour levels). Furthermore, all the modes are located at the correct 
frequencies (within Au ~ 0.125). For At = 0.5, one has cjat = 27r thus normally 
allowing for modes n = and n = 1 only. In the map one can recognize the 
n = 0, 1,2 modes. None of them is located at the correct frequency, it is also 
surprising to find the n = 2 mode whose eigenfrequency is above the Nyquist 
frequency. In this case it is clear that the time step is too large to accurately 
resolve the sound waves. To analyze the gravity modes we use the method devel- 
oped in DB: we project the velocity field on the anelastic eigenvectors. With this 
method we obtain the individual mode amplitudes q„ which are the coefficients 
in the eigenvector expansion. With such a decomposition it is straightforward to 
obtain the mode frequency. We have checked (not shown here) that for At = 0.5 
(run 3), the gravity mode frequencies perfectly match the analytical predictions. 

This test highlights the importance of the choice of the time step to cor- 
rectly describe the physical process of interest. In the present case, sound wave 
amplitudes are much lower than the amplitudes of the gravity modes (the ini- 
tial setup is almost in hydrostatic equilibrium) and therefore they do not play 
an important role in this problem, one can then safely use a larger time step 
(e.g. At = 0.5) to study the internal gravity modes, as shown by the excellent 
agreement with predicted frequencies. Obviously one could not do that if the 
initial setup was intended to be far from hydrostatic equilibrium as sound waves 
would then be at least as important as internal gravity modes. 



Figure 2. Vertical profiles of the horizontally averaged fluxes (full: total 
flux, dashed: convective flux, dotted: radiative flux) for m = 1, 0, —0.9 (from 
the left to the right). The expected ratio of i^conv/^tot are 20%, 60% and 
96%. 



3.2. Test 2: Rayleigh-Benard convection 

Since one of the first planed astrophysical applications of our implicit code is 
stellar convection, the simple problem of Rayleigh-Benard convection is a good 
test. The initial setup and parameters are taken from Brandenburg et al. (2005). 
We consider a stratified atmosphere in a cartesian domain [—0.5,0.5] x [0,0.5]. 
The initial atmosphere is a polytropic model characterized by a polytropic index 
m so that initially p oc T™. We consider a perfect gas with 7 = 5/3. For 
m < 1/(7 — 1) = 3/2 the temperature stratification is super-adiabatic and 
the atmosphere is dynamically unstable. If in addition the Rayleigh number is 
above its critical value, the atmosphere will be also convectively unstable. In 
Brandenburg et al. (2005) it is shown that for a given value of m, one can 
deduce the expected ratio of the convective flux Fconv to the total flux Ftot • For 
example, for m =1, and —0.9, this ratio is respectively 20%, 60% and 96%. 
The purpose of this test is to check if we can reproduce these values. 

We use a constant viscosity v = 5x 10~^ and a constant thermal conductiv- 
ity k = 10~^. Periodic boundary conditions are used in the horizontal direction. 
Non penetrative conditions are used at the top and at the bottom where in 
addition the total flux Ftot is imposed. The Rayleigh number is computed as 
Ra = jr^^ where x = k/{poCy) is a mean value of the thermal diffusivity 
(To and po are mean values of the temperature and density). All simulations 
are done on a 50^ grid. All run parameters are summarized in Table [2j In each 
run the system converges to a steady state showing two convective cells. This 
is not surprising since the aspect ratio of our domain is chosen so that the most 
unstable mode almost exactly fits in the box. This steady state is reached after 
about 100 units of time. In practice the time step converges at some value shown 



Table 2. Run parameters for the Rayleigh-Benard test 



Run 


m 


Ra 


Resolution 


max. time step 


max. CFL 


1 


1 


5 X 10^ 


50^ 


0.04 


1 


2 





10^ 


502 


0.5 


12.5 


3 


-0.9 


10^ 


502 


0.5 


37.5 
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Figure 3. Snapshot showing the convective eddies. The axis are normalized 
by the stellar radius R^,. The color codes the logarithm of the density (light 
to dark color means high to low density). 

in Table El The m = —0.9 case has a smaller time step, which is not surprising 
since it is the most unstable case investigated here. 

We define the following horizontally averaged fluxes -Frad =< >> 

2 _ 

^enth =< pUzCp{T- < T >) >, Fkin =< pUz\ > and Fvisc =< -u.f >, where 
< . > stands for horizontal averaging. In the steady state, their sum should be 
equal to -Ftot- In practice we find that this is the case within 0.5%. The kinetic 
energy fiux i<kin and the viscous flux -Fvisc are found to transport less than ~ 1% 
of the total flux so that Fconv ~ -^enth • Figure [2] displays the vertical profiles of 
these fluxes and shows that the expected ratios are recovered. Note that in the 
last case, almost all the flux is transported by convection. 

4. Preliminary results under realistic stellar interior conditions 

We consider a hot, almost fully radiative star model to test the ability of the 
code to describe convection in stellar conditions. The initial stellar model used 
for the hydrodynamic calculations is produced by a ID stellar structure code 
and has T^a = 7500 K, M = 5Mq, log{L/LQ) = 3 and log{g) = 2.6. The mixing 
length parameter (hereafter MLT) is taken as 1.7. This model lies outside the 
instability strip of Cepheids so it is stable against pulsations. This model has a 
small convective region near the surface, in the H/He ionization region around 
T ~ 10^ K. Our initial model covers 60% of the stellar radius i?^ ~ 19Rq. 

We consider a numerical domain (r, 9) G [0.4i?^, i?^] x [vr/S, 2tt/3]. We use a 
restricted angular domain to increase the angular resolution without sacrifying 
too much CPU time. The grid size is 150^. We use reflective boundary conditions 
in the 6 direction and non-penetrative conditions at the top and at the bottom of 
the envelope. An energy flux consistent with the stellar luminosity is imposed at 
the bottom. The boundary conditions for the energy flux at the outer boundary 
are not straightforward and different solutions are possible. As a preliminary 
solution, we impose that the last layer of cells radiates toward the exterior a 
flux equal to uT^, where T is the temperature of the cells. This is justified by 
the fact that the outer boundary conditions of the initial ID model relies on the 
Eddington approximation and is defined at the photosphere, at an optical depth 
TRoss = 2/3 where T = Teg. 

The initial ID model does not exactly correspond to an equilibrium state 
of the hydrodynamic code. Consequently, the model undergoes an initial phase 
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Table 3. Comparison with MLT. 





Simulation 


MLT 


Convection Mach number Vconv/cs 


0.02 


0.05 


Eddy turnover time-scale 


lO^s - lO'^ s 


10^ s 


Eddy size 


~ 5Hp 


1.7 Hp 


-^conv/-^tot 


3 X 10"^ 


2 X 10^3 



of relaxation, for which we find more convenient to use a dissipative, more ro- 
bust, 1st order scheme (Euler backward + donor cell). This initial relaxation 
phase lasts for about 100 days during which the model reaches hydrostatic and 
thermal equilibrium, which corresponds to the Kelvin-Helmholtz time-scale of 
the envelope. Convective eddies then develop near the surface (see Fig. [3|). At 
this stage, we switch to a second order scheme. The eddies remain throughout 
the whole simulation, which lasts over 700 days. The mean time step is ~ 10^ s 
which corresponds to a CFL number of 400. 

The eddies pattern show a non-steady behavior, but always remain confined 
near the top boundary. In table [3l it is shown that our results are consistent 
with MLT. It is important to note that the mean time step of the simulation 
is lower by a factor of ten than the maximum eddies turnover time-scale. This 
tells us that the latter quantity is certainly relevant to determine the time step 
in order to accurately describe convection. 

5. Conclusion 

We have presented a new multi-D hydrodynamic implicit code. As shown in 
section 13. H high CFL simulations are efficient to study physical processes pro- 
ceeding on time-scale larger than the CFL limit. One should always consider 
numerical results with care, as an inappropriate time step can give significantly 
wrong results. We have shown that the code is able to handle strong convection, 
at least in a simple setup. Our first results on a realistic (weakly convective) 
stellar model are encouraging and provide very promising perspectives for the 
first multi-D simulations of convection in pulsating star, a challenge in the field 
of asteroseismology. 
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